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This  report  presents  the  motivation  and  mathematical  background 
for  the  parametric  terrain  model  used  In  the  STAR  Brigade  level  combined 
arms  combat  simulation.  Computer  subroutines  for  terrain  elevation  and 
line  of  sight  computations  are  presented  and  explained  In  detail  along 
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I.  INTRODUCTION 

A.  Origins  of  the  Parametric  Terrain  Model 

The  STAR  (Simulation  of  Tactical  Alternative  Responses)  ground-air 
combat  model  Is  an  extensive  computer  simulation  program  developed  at  the  Naval 
Postgraduate  School  during  1978-1979.  STAR  is  written  in  the  SIMSCRIPT  II. 5 

simulation  language.  This  report  documents  the  battlefield  terrain  representa- 
tion which  Is  used  In  the  current  Brigade-level  version  of  the  model.  The  idea 
for  this  terrain  representation— called  parametric  terrain--was  originally  pro- 
posed by  Maj.  Chris  Needels  in  his  1976  Master  of  Science  thesis  at  NPS  [ 1 ]. 

The  basic  function  which  any  terrain  representation  must  provide 
for  a high-resolution  combat  simulation  Is,  "for  any  X,  Y map  coordinates  on 
the  battlefield,  compute  the  elevation  Z of  the  terrain."  This  elevation  Z 
Is  generally  called  the  macro  terrain.  Macro  terrain  provides  a somewhat  smoothed 
replica  of  an  actual  battlefield  in  that  very  small  features  (e.g.  a 1 meter 
boulder)  are  not  represented,  but  major  features  (such  as  a 200  m 'll)  are 
represented. 

B.  Approaches  to  Terrain  Modelling 

Most  current  high  resolution  combat  simulations  use  a macro  terrain 
representation  known  as  "digitized  terrain".  Essentially  digitized  terrain  in- 
volves storing  a (usually  large)  table  of  elevations,  Z,  for  a grid  of  X,  Y 
coordinates  covering  the  battlefield.  The  process  of  determining  Z for  a given 
X,Y  then  reduces  to  one  or  more  table  look-ups,  possibly  followed  by  an  inter- 
polation to  smooth  the  terrain  between  grid  points.  Digitized  terrain  enjoys 
simplicity  and  speed  of  computation,  but  requires  extensive  computer  storage  de- 
voted to  the  elevation  table.  This  has  tended  to  restrict  the  size  of  the 
battlefield  for  models  using  digitized  terrain  to  about  10  x 10  km. 

The  parametric  terrain  model  used  in  STAR  essentially  involves  storing 
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a (rather  complicated)  function  f(X,Y).  The  process  of  determining  Z for 
a given  X,Y  then  reduces  to  computing  the  function  Z * f(X,Y).  Parametric 
terrain  has  the  advantage  that  the  function  f can  be  represented  using  only 
a modest  amount  of  computer  storage,  thus  allowing  simulations  to  be  run  on 
much  larger  battlefields.  In  addition,  the  parametric  terrain  Is  inherently 
continuous,  so  no  Interpolation  Is  required  for  smoothing,  and  slopes  can  be 
computed  directly  with  no  difference  approximations  required.  The  associated 
computations,  however,  may  be  quite  complex,  and  great  care  must  be  taken  to 
achieve  computational  efficiency. 

Since  STAR  Is  a Brigade  level  combined-arms  simulation,  with  long  range 
systems  (air,  artillery,  resupply)  represented  on  the  battlefield,  it  is  im- 
portant to  have  a large  terrain  area  available.  Thus  STAR  incorporates  a 
parametric  macro  terrain  model. 

Discussions  with  experienced  modelers  have  Indicated  that  attempts 
have  been  made  In  the  past  to  develop  terrain  models  which  share  some  of  the 
features  of  the  STAR  parametric  terrain  model.  These  efforts  do  not  seem  to 
* have  been  successfully  Incorporated  Into  combat  models,  but  since  documentation 
Is  generally  not  available,  the  reasons  for  this  are  not  clear. 

C.  Overview  of  the  Report 

Section  II  of  this  report  will  develop  the  parametric  terrain  rep- 
resentation and  will  Indicate  how  It  can  be  used  to  emulate  actual  terrain  areas. 
In  Section  III  we  detail  the  process  of  computing  elevations,  including  the 
STAR  routines  ELEV  and  ELEVG.  A terrain  preprocessor  program  HILL. LIST  for 
Increasing  the  efficiency  of  computations  is  also  discussed.  Section  IV  dis- 
cusses augmenting  the  terrain  representation  to  Include  forested  areas. 

Given  a terrain  representation,  one  vital  combat  function  which  must  be 


computed  Is  the  existence  or  absence  of  line-of-sight  (LOS).  Section  V 
discusses  In  general  terms  the  STAR  LOS  procedure.  Section  VI  goes  into  the 
mathematics  of  the  LOS  computations  in  substantial  detail,  and  analyzes  the 
computer  code. 

Section  VII  presents  the  computer  code  for  the  terrain  and  forest 
preprocessor  programs  and  for  the  data  read-in  program  for  the  terrain. 

Finally,  in  Section  VIII  we  discuss  some  areas  which  seem  to  have 
potential  for  further  research  in  parametric  terrain  modelling. 

II . The  STAR  Parametric  Macro  Terrain  Model 
A,  General  Form  of  the  Terrain  Model 

The  parametric  terrain  model  proposed  by  Needels  in  [1]  represents 
terrain  by  modelling  individual  hill  masses.  The  overall  terrain  is  then 
obtained  by  superpositioning  the  individual  hills.  Mathematically,  if  fj(X,Y) 
is  a function  giving  the  elevation  of  the  Ith  hill  mass  at  any  X,Y,  then  the 
overall  terrain  elevation  at  X,Y  is  obtained  as  the  pointwise  maximum  over 
all  the  hill  masses, 

Z = f(X,Y)  * maximum  fj(X,Y)  (1) 

1=1,2,. ...NHILLS 

where  NHILLS  is  the  total  number  of  hill  masses  on  the  battlefield. 

A schematic  cross-section  view  of  several  hill  masses  along  with 
the  resulting  terrain  elevation  is  shown  in  Figure  1. 

In  the  Needels  model  each  individual  hill  mass  is  represented 

\ 

mathematically  as  a scaled  bivariate  normal  probability  density  function. 

This  gives  a characteristic  bell-shaped  hill  mass  cross-section  as  shown  in 

Figure  1,  and  elliptical  contours  for  each  individual  hill  mass.  By  varying 

the  bivariate  normal  parameters,  a wide  variety  of  different  hill  locations, 

sizes,  and  shapes  can  be  modelled.  By  superpositioning  several  hill  masses, 

the  contour  map  can  be  fit  to  real  map  contours  remarkably  closely. 
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The  STAR  parametric  terrain  model  uses  a generalization  of  the  hill 
mass  functions  of  the  Needels  model.  The  representation  of  these  functions  has 
been  changed  significantly  to  make  them  easier  to  fit  to  real  terrain  maps  and 
to  make  computations  more  efficient.  In  addition,  forest  features  have  been 
added,  and  the  rather  complex  llne-of-slght  computations  have  been  developed. 

B.  Parameterization  for  Terrain  Fitting 

In  order  to  emulate  a piece  of  real  terrain  with  the  STAR  parame- 
tric terrain  model.  It  Is  necessary  to  fit  the  model's  hill  mass  functions 
f j ( X , Y ) to  a contour  map  of  the  terrain  to  be  modeled.  The  fitting  process  is 
currently  done  by  hand,  so  It  is  essential  to  describe  the  hill-mass  functions 
using  a parameterization  that  makes  geometric  sense  and  can  easily  be  related 
to  a contour  map.  In  this  section  we  develop  the  "fitting-parameterization"  of 
the  hill  mass  functions  fj(X,Y).  This  parameterization  is  not  particularly 
well  suited  for  efficient  computation,  so  we  will  convert  to  a computing- 
parameterization  in  Section  III. 

Each  hill  mass  function  fj(X,Y)  is  given  by  an  exponential 

function 

f j ( X , Y ) = PEAK. H( I ) + HT.H(I)*  [exp  (Q^X.Y))  - 1]  (2) 

where  PEAK.H(I)  and  HT.H(I)  are  constant  parameters  for  hill  I to  be  described 
shortly,  and  where  Qj ( X* Y)  is  a quadratic  function  of  X,Y  which  has  several 
other  parameters.  Qj(X,Y)  is  a negative  definite  quadratic,  so  the  exp(Qj(X,Y)) 
term  yields  hills  having  the  characteristic  bell-shaped  cross-section  of  a 
normal  probability  density  function,  and  having  contour  lines  which  are  elliptical. 
The  parameters  control  the  location,  orientation,  size,  and  shape  of  the  ellip- 
tical contours. 

Elevations  in  the  model  are  measured  from  0 meters  = sea-level,  and 
the  terrain  as  a whole  has  a minimum  elevation  level  above  sea  level  denoted  BASE. 
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The  BASE  value  Is  Included  as  a tern  in  the  maximization  of  equation  (1) 

Z * max  (BASE,  max  fT(X,fj>  (3) 

I 1 

Several  of  the  hill  parameters  relate  to  hill  height,  and  are  best 
visualized  on  a cross-section  diagram  (see  Figure  2.). 

HT.H(I):  the  maximum  height  of  the  "normal" 

curve  describing  this  hill  mass. 

PEAK.H(I):  the  elevation  of  the  hilltop  measured 
from  zero  = sea-level. 

Since  HT.H(I)  may  be  greater  than  or  less  than  PEAK.H(I),  the  bottom  of  the 
hill  which  occurs  at  PEAK.H(I)  - HT.H(I)  need  not  be  at  sea-level.  However, 
since  normal  curves  have  tails  approaching  the  botom  level  which  extend  to  + ® 
we  require  that  these  tails  be  below  the  terrain  BASE  value.  Thus  a constraint 
on  the  terrain  parameterization  is 

PEAK . H ( I ) - HT.H(I)  < BASE  (4) 

For  some  hill  masses,  the  tails  of  the  hill  are  annoying  as  they  pre- 
vent the  terrain  from  dropping  off  quickly  enough  to  model  a cliff  side  or  a 
steep  riverbank.  Thus  we  include  the  parameter 

CUT.H(I):  A vertical  distance  measured  down  from 
the  peak  beyond  which  this  hill  mass  is 
no  longer  considered  in  the  computations. 

Thus  if  CUT.H(I)  is  70m.,  at  most  seven  contours  of  hill  I will  appear 
on  a 10m.  contour  map.  The  CUT.H  parameter  must  be  used  with  care,  since  it 
may  introduce  unwanted  discontinuities  in  the  terrain--vertical  cliffs.  Figure 
3 shows  an  example  in  which  this  parameter  is  properly  used.  Here  the  gradual 
slope  of  hill  2 has  been  terminated  using  CUT . H (2 ) and  hill  3 has  been  added 
to  give  a more  abrupt  slope  to  the  terrain  at  this  point.  Without  hill  3 the 


* 
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resulting  terrain  would  have  a bad  discontinuity  as  in  Figure  4.  For  *he  vast 
majority  of  the  hills  not  requiring  the  CUT.H(I)  parameter  it  may  be  safely 
defaulted  to  HT.H(I).  The  CUT.H  parameter  also  improves  computational  effi- 
ciency by  restricting  the  region  in  which  we  need  to  worry  about  hill  I.  The 
Terrain  Preprocessor  program  HILL. LIST  to  be  described  in  Section  III  will 
compute  an  Improved  CUT.H  value  for  each  hill  (but  will  be  careful  not  to 
introduce  discontinuities). 

The  remaining  fitting  parameters  are  most  easily  understood  in  reference 
to  the  elliptical  contour  map  of  a hill  mass.  Figure  5 shows  the  top  view 
contour  plot  of  hill  mass  I,  with  only  one  contour  line  shown  50m.  down  from 
the  peak, to  avoid  cluttering  the  diagram.  The  following  parameters  are  indicated: 

XC.H(I):  the  X map-coordinate  of  the  center  location  of  hill  I 

YC.H(I):  the  Y map-coordinate  of  the  center  location  of  hill  I 

ANG.H(I):  the  orientation  angle  of  the  ellipse  measured  in  degrees 

counter-clockwise  from  EAST,  to  the  major  axis. 

(in  Fig.  5,  ANG-H  = 30°) 

ECC.H(I):  the  eccentricity,  defined  as  the  ratio  of  major  axis 

length  to  minor  axis  length.  (ECC.H  _>  1) 

(in  Fig.  5,  ECC  = 2.0) 

SPRD.H(I):  spread  is  a measure  of  hill  size  which  is  defined  as 

the  distance  in  meters  measured  along  the  major  axis 
from  hill  center  to  the  contour  line  which  is  50  meters 
down  from  the  peak. 

One  consequence  of  this  somewhat  arbitrary  definition  of  the  spread  parameter 
is  that  all  hills  must  have  a HT.H(I)  greater  than  50m. 
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. Fitting  Hills  to  Real  Terrain 


Given  the  set  of  8 hill  parameters,  (XC.H,  YC.H,  PEAK.H,  ANG.H, 

ECC.H,  SPRD.H,  HT.H,  CUT.H),  our  goal  Is  to  emulate  a piece  of  real  terrain 
by  setting  the  parameters  to  mimic  a real  contour  map.  The  process  is  an 
Iterative  manual  fit- & -evaluate  cycle  which  requires  some  practice  and  the 
development  of  an  Intuition  relating  the  hill  parameters  to  the  geometric 
notions  of  contour  shape. 

Preliminary  analysis  of  a terrain  area  indicates  the  number, 
location,  and  peak  of  the  primary  hill  masses,  and  a guess  at  the  number  of 
parametric  hills  needed  to  represent  each  real  hill  mass.  Then  the  parameters 
for  each  hill  are  set,  and  a computer  generated  contour  map  is  drawn  and 
compared  to  the  original.  Discrepancies  are  noted,  and  parameters  adjusted  to 
improve  the  fit.  After  a few  weeks  of  practice,  our  experience  shows  that  an 
Individual  can  fit  a 10  x 10  km  map  section  In  about  a week. 

The  resulting  parametric  terrain  can  be  made  remarkably  similar  to 
the  original  terrain  as  Is  shown  in  the  comparative  contour  maps  of  Figures  6-11. 


As  a rough  Indication  of  the  computer  storage  economies  which  are 
possible  with  this  terrain  parameterization,  if  a 10  km  x 10km  battlefield 
can  be  represented  using  NHILLS  = 100  hills,  then  800  parameters  must  be  stored. 
By  contrast,  a digitized  terrain  model  working  with  a 100m  grid  size  would  have 
to  store  slightly  more  than  10,000  numbers.  This  represents  a storage  savings 
factor  of  about  12.  A more  careful  comparison  would  have  to  consider  the  exact 
number  of  hills  required  as  well  as  possibilities  for  storage  packing  In  both 
model s . 
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III.  Elevation  Computations 
A.  Parameterization  for  Computing 

The  fitting  parameters  discussed  in  Section  II  were  designed  for 
their  geometric  connotations  to  aid  in  the  terrain  fitting  process.  For  com- 
puter calculations  we  are  Interested  in  computational  efficiency  rather  than 
geometric  appeal,  so  a second  parameter  set  has  been  developed  for  the  STAR 
terrain  model.  This  set  of  computing-parameters  is  derived  algebraically  from 
the  fitting  parameters  and  yields  exactly  the  same  terrain  elevations. 

The  basic  elevation  formula  for  hill  I is,  as  in  equation  (2), 


fj(X.Y)  * PEAK.H(I)  - HT.H(I)  + HT.H(I)*exp  [Qj(X.Y)]  (5) 

where  PEAK.H(I)  - HT.H(I)  gives  the  elevation  of  the  bottom  of  the  hill, 
and  HT.H(I)  *exp  [Qj(X,Y)]  adds  the  hill  mound  to  this  bottom  level  as  a 
function  of  X and  Y . 

For  the  computing  parameterization,  we  define 


XS  = X - XC.H(I) 
YS  = Y - YC.H(I) 


(6) 


giving  the  coordinate  distances  from  X,Y  to  the  hill  center  location.  Then 
the  quadratic  function  Q j ( X , Y ) is  given  as 

Qj(X,Y)  * PXX.H(I)*XS2  + PYY.H(I)*YS2  + PXY.H(I)*XS*YS  (7) 

where  PXX.H(I),  PYY.H(I),  and  PXY.H(I)  are  computing-parameters  defined 
in  terms  of  the  fitting  parameters  as  follows: 

Letting  SANG  « sin  (ANG.H(I)) 

CANG  « cos  (ANG.H(I)) 

A * 

B - A * ECC ( I V 


In  (HT.H( I )/(HT.H( I ) - 50.)) 


(8) 
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then  PXX.H(I)  =■  - (A*CANG2  ♦ B *SANG2)/SPRD.H(I  )2 

PYY.H(I)  = - (A*SANG2  + B *CANG2 ) /SPRD . H ( I )2  (9) 

PXY.H(I)  - (2*CANG  *SANG  * (B-A) )/SPRD. H( I )2 
Finally^  the  CUT.H(I)  parameter,  which  has  been  ignored  in  all  of  the  above 
formulas  -is  related  to  Qj(X,Y)  by  defining  the  computing  parameter 

(In  ((HT.H(I)  - CUT.H( I ) )/HT.H( I ) ) 

if  CUT.H(I)  < HT.H(I)  - (10) 

- » otherwise 

and  redefining  the  basic  hill  equation  to  be 


f,(X,Y) 


PEAK.H(I)  - HT.H(I)  + HT.H(I)*exp  [Qj(X.Y)] 

if  Qj (X ,Y)  > CRIT.H(I)  (11) 

BASE  otherwise. 


The  derivation  of  these  formulas  is  rather  tedious,  but  certainly 
not  profound,  consisting  primarily  of  a coordinate  rotation,  the  equation 
of  an  ellipse,  and  the  definitions  of  SPRD.H  and  CUT.H.  The  reader  who  wishes 
to  ful^y  understand  the  mathematics  of  the  terrain  model  should  work  through 
these  derivations. 

B.  A Primitive  ELEV  Routine 

Given  the  above  computing  parameters  (XC.H,  YC.H,  PEAK.H,  HT.H , 
PXX.H,  PYY.H,  PXY.H,  CRIT.H)  for  each  of  a number  NHILLS  of  hills  on  our 
battlefield,  and  given  an  overall  terrain  BASE,  a primitive  computer  subroutine 
can  easily  be  written  to  compute  the  elevation  Z for  any  point  X,Y  on  the 
battlefield.  Such  a program  is  given  here  in  the  SIMSCRIPT  II. 5 computer 
simulation  language,  assuming  the  above  parameters  have  already  been  stored  in 
appropriate  global  arrays  and  variables. 


ROUTINE  FOR  PRIMITIVE. ELEV  GIVEN  X AND  Y YIELDING  Z 

NORMALLY  MODE  IS  REAL 

DEFINE  I AS  AN  INTEGER  VARIABLE 

LET  Z = BASE 

FOR  1=1  TO  NHILLS  DO 

LET  XS  = X - XC.H(I)  LET  YS  = Y - YC.H(I) 

LET  QI  = PXX.H(I)  *XS**2  + PYY.H(I)  *YS**2  + PXY.H(I)*XS*YS 
IF  QI  LT  CRIT.H(I)  CYCLE 

ELSE  LET  FI  = PEAK.H(I)  + HT. H ( I )*( EXP . F (QI ) - 1 . ) 

IF  FI  GT  Z LET  Z = FI  ALWAYS 
LOOP  RETURN  END 

Given  the  formulas  in  Section  III. A,  this  routine  should  not  require 
further  documentation.  It  is  the  protopype  for  the  ELEV  routine  actually 
used  in  STAR. 

C.  Refinements  for  Speeding  the  Calculations 

The  PRIMITIVE. ELEV  routine,  while  it  will  work  for  any  X,Y,  is 
not  particularly  efficient.  Its  primary  weakness  is  the  need  to  loop  over  all 
NHILLS  hills  even  though  only  a few  hills  are  close  enough  to  a given  X,Y  to 
have  any  chance  of  influencing  the  macro-terrain  elevation  Z at  that  point. 
Substantial  economies  in  computing  time  can  be  realized  by  using  a terrain 
preprocessor  program  to  develop  lists  of  the  hills  which  are  relevant  to  par- 
ticular areas  of  the  battlefield  and  to  improve  CUT.H  values  for  all  hills. 

The  terrain  preprocessor  program,  called  HILL. LIST,  divides  the  battlefield  into 
grid  squares  of  size  6SIZE  (our  experience  indicates  that  GSIZE  in  the  range 
1 to  3 km  yields  an  effective  grid).  Each  grid  square  is  referenced  by  grid 
subscripts  IX  and  IY.  An  array  LIST.H  (IX,IY,L)  is  created  which  contains, 
for  each  IX,IY  (i.e.  for  each  grid  square)  a list  of  the  hills  I which 
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contribute  to  the  terrain  in  the  grid  square. 

Then  given  X,Y  we  can  compute  the  elevation  Z by  first  computing 
the  appropriate  grid  square  IX, IY  and  then  executing  a routine  similar  to 
PRIMITIVE. ELEV  which  only  loops  over  the  hills  in  LIST. H( IX , IY ,L ) for  this 
grid  square.  Dramatic  computational  efficiencies  result  with  only  a modest 
storage  increase  for  the  list.  Typical  results  from  our  experience  are  that 
for  6S1ZE  = 1 km. , on  the  average,  4 or  5 hills  are  relevant  to  each  grid 
square.  Clearly  this  depends  on  the  complexity  of  the  terrain  and  the  number 
of  hills  used  to  model  it. 

A further  savings  can  be  obtained  by  computing  in  HILL. LIST  the 
elevation  E of  the  lowest  terrain  point  actually  represented  by  each  hill 
and  using  E with  PEAK.H  to  get  the  smallest  CUT.H  value  possible  for  the 
hill.  Then  the  QI  vs  CRIT.H(I)  test  in  ELEV  will  cycle  more  frequently  hence 
avoiding  the  expensive  exponential  computation. 

The  line-by-line  details  of  HILL. LIST  are  documented  in  Secion  VII. 

For  this  section  it  suffices  to  define  the  global  LIST.H  (IX,IY,L)  array  which 
HILL. LIST  creates.  LIST.H  has  subscripts 

IX:  grid  square  subscript  in  X direction 

IX  = 1 , . . . .NGRIDX 

IY:  grid  square  subscript  in  Y direction 

IY  = 1.....NSRIDY 

L:  subscript  for  last  coordinate  of  LIST.H  . 

For  a given  IX, IY, 

LIST. H( IX, IY , 1 ) = BASE  value  for  the  grid  square 

L 1ST. H ( IX ,{Y,2) , LIST. H(IX,IY,3) ,ETC.  = hill  numbers  for  the  hills 
which  influence  terrain  in  this  gridsquare. 


LIST.H  Is  represented  In  SIMSCRIPT  as  a ragged  array,  so  the  number  of  hills 
to  be  scanned  In  grid  square  IX, IY  Is  DIM.F(LIST.H(IX,1Y,*))-1  . 

The  resulting  streamlined  ELEV  routine  Is  presented  In  Listing  1. 

Local  variables 

I hill  number 
ix7 

jy  I square  Index  In  X and  Y directions 

KOUNT  number  of  hills  in  a grid  square  + 1 
L loop  index 

Y ^ Input  X and  Y coordinates 
Z resulting  elevation 

XS,YS,QI,FI  intermediate  computations  as  In  PRIMITIVE. ELEV 
Global  variables 
X LO  BDRY^ 

Y^LO.BDRY  i map  c00rc,1nates  southwest  corner  of  battlefield  in  meters 
GSIZE  grid  square  size  In  meters 

NGRIDX  range  of  IX  - number  of  grid  squares  In  X direction 

NGRIDY  range  of  IY  - number  of  grid  squares  in  Y direction 

LIST.H  list  of  hill  numbers  for  each  grid,  also  BASE  value 

DUM.I  one-dimensional  integer  dummy  array  to  simplify 

coordinate  computation  in  accessing  LIST.H 
XC.H,YC.H,PEAK.H,HT.H,PXX.H,PYY.H,PXY.H,CRIT.H  hill  parameter  arrays. 
Routines  called 
None 

Events  scheduled 
None 
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Lines  11-12 

Line  13 
Lines  14-20 


Lines  5-10  compute  the  grldsquare  Indices  and  make  sure  that  the  grid- 
square  Is  on  the  battlefield, 
define  DUM.I  to  be  the  appropriate  column  of  LIST.H 
and  determine  the  length  of  that  column 
sets  Z to  the  BASE  value 

loop  over  the  hills  I (for  this  grid  square  only) 
and  Increase  Z if  the  computed  hill  elevation  FI  is 
greater  than  the  max  so  far. 

Note  that  the  origin  for  the  grid  square  reference  system  is  at 
X.LO.BDRY,  Y.LO.BDRY  (the  southwest  corner  of  the  battlefield). 

D.  Terrain  Slope 

In  addition  to  terrain  elevation,  the  slope  of  the  terrain  at  any 
point  X,Y  is  directly  available.  Differentiating  equation  (11)  gives  the 
gradient  components 

HT.H(  I )*exp[QT (X  ,Y)]*  8QI(X,Y) 


GX  = 3fj(X,Y) 

aT — 


where 


3Qj(X,Y) 

ax 


and 


GY  = »VX*Y) 
aY 


where 


3 QjtX.Y) 
aY 


aX 

if  Qj(X.Y)  > CRIT.H(I) 
0 otherwise 

2*PXX.H( I )*XS  + PXY . H( I )*YS 


HT.H(  I )*exp[Qj  ( X » Y ) ] * 


If  Qj(X.Y)  > CRIT.H(I) 
otherwise 


2*PYY.H(I)*YS  + PXY.H(I )*YS  . 


(12) 


(13) 


(14) 


(15) 


Then  If  we  move  away  from  X,Y  In  direction  AX,  AY  the  slope  Is  given 
by  the  directional  derivative  / 

D ■ (GX  * AX  + GY  * AY)/yAX2  + AY2  . (16) 

An  elevation  routine  called  ELEVG  which  computes  the  gradient  components 
as  well  as  Z Is  available  and  is  used  by  the  STAR  movement  routines.  A 
separate  routine  Is  used  because  the  line  of  sight  routine  which  calls  ELEV 
does  not  need  the  gradient  information  and  should  not  pay  the  added  computational 
price.  The  computational  sequence  In  ELEVG  Is  slightly  different  to  capitalize 
on  common  computations  in  Z,  GX,  and  GY.  The  resulting  code  is  given  in  Listing  2. 

Because  of  the  similarity  to  ELEV  we  will  not  provide  a line-by-line  explan- 


ation of  ELEVG. 


IV.  Forest  Model  line 


A.  Brief  Description  of  Forest  Methodology 

In  addition  to  the  macro-terrain,  another  factor  which  crucially 
influences  line  of  sight  computations  In  the  STAR  model  (as  in  the  real  world) 

Is  the  presence  of  forested  areas.  Although  we  will  refer  to  "forests" 
throughout  this  discussion,  the  modeling  tools  developed  here  could  just  as  well 
be  applied  to  man-made  features,  such  as  towns,  which  provide  cover  and  thus 
interrupt  LOS. 

The  same  basic  ideas  can  also  be  expanded  to  model  smoke  clouds  on  the 
battlefield.  The  resulting  clouds  must  be  able  to  move,  expand,  and  disperse 
(which  forests  do  not  do),  and  are  not  totally  opaque  to  all  sensor  systems. 
Development  of  a smoke  module  for  STAR  is  currently  underway,  but  It  will  not 
be  documented  here. 

Forests  or  other  cover  features  In  the  STAR  model  are  represented  by 
elliptical  areas  on  the  ground.  Each  such  cover  ellipse  has  a tree  height 
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Listing  2 


associated  with  It,  and  the  forest  Is  thus  an  elliptical  "cylinder"  with 
that  fixed  height  above  the  macro  terrain.  Forests  with  non-elliptical  shapes 
and  non-constant  heights  can  be  approximated  by  combining  several  possibly 
overlapping  ellipses.  The  tree  height  at  a given  point  X,Y  Is  the  maximum 
tree  height  for  all  the  forest  ellipses  containing  the  point  X,Y. 

B.  Parameterization  for  Forest  Fitting 

As  Is  the  case  for  macro  terrain  hill  fitting,  a geometrically 
motivated  parameterization  of  the  cover  ellipses  has  been  developed  to  aid 
In  the  fitting  of  the  model  to  terrain  maps.  The  parameters  are  listed  below 
and  Illustrated  In  Figure  12. 

XC.E(I):  the  X coordinate  of  the  ellipse  center 

YC.E(I) : the  Y coordinate  of  the  ellipse  center 

ANG.E(I):  the  orientation  angle  measured  counterclockwise 

in  degrees  from  east  to  the  major  axis  of  the  ellipse 
AMAJ.E(I):  length  In  meters  of  the  semi-major  axis  of  the  ellipse 
AMIN.E(I):  length  of  semi-minor  axis  of  the  ellipse 
HT.E(I):  height  of  the  trees  In  ellipse  I above  the  macro  terrain 

elevation. 

C.  Parameterization  for  Computing 

Given  an  arbitrary  point  X,Y  on  the  battlefield  we  want  to  de- 
termine If  this  point  Is  Inside  a forest  feature,  and.  If  so,  how  high  the 
trees  are.  Since  this  computation  will  be  performed  frequently  in  the  LOS 
procedure,  it  Is  important  to  have  an  efficient  computational  procedure.  Thus 
the  ellipses, like  the  hills,  have  a transformed  parameter  set  for  computing. 


The  boundary  of  the  Itn  forest  ellipse  Is  represented  by  the  quadratic 


equation 
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tree  height  for  all  the  forest  ellipses  containing  the  point  X,Y. 

B.  Parameterization  for  Forest  Fitting 

As  Is  the  case  for  macro  terrain  hill  fitting,  a geometrically 
motivated  parameterization  of  the  cover  ellipses  has  been  developed  to  aid 
in  the  fitting  of  the  model  to  terrain  maps.  The  parameters  are  listed  below 
and  Illustrated  In  Figure  12. 

XC.E(I):  the  X coordinate  of  the  ellipse  center 

YC.E(I):  the  Y coordinate  of  the  ellipse  center 

ANG.E(I):  the  orientation  angle  measured  counterclockwise 
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C.  Parameterization  for  Computlnc 


Given  an  arbitrary  point  X,Y  on  the  battlefield  we  want  to  de- 
termine If  this  point  is  Inside  a forest  feature,  and.  If  so,  how  high  the 
trees  are.  Since  this  computation  will  be  performed  frequently  In  the  LOS 


procedure.  It  Is  Important  to  have  an  efficient  computational  procedure.  Thus 
the  ellipses, like  the  hills,  have  a transformed  parameter  set  for  computing. 

The  boundary  of  the  Ith  forest  ellipse  Is  represented  by  the  quadratic 


equation 


(17) 


PXX.E(I)*(X-XC.E(I))2  ♦ PYY.E(I)*(Y-YC.E(I))2 
♦ PXY.E(I)*(X-XC.E(I))*(Y-YC.E(I))  - 1.0 

Then  an  arbitrary  X,Y  Is  inside  the  ellipse  If  and  only  If  the  left  hand 
side  of  equation  (17)  Is  less  than  1 . 

In  terms  of  the  fitting  parameters  the  quadratic  equation  coefficients 
are  given  as: 

SANG  - sin  (ANG.E(I)V 
CANG  ■ cos  (ANG.E(I))  . 

PXX. E(I)  - (CANG/AMAJ.EO))2  ♦ ( SANG/AMIN. E(I))2  (18) 

PYY.E(I)  - (SANG/AMAJ.E(I))2  ♦ (CANG/AMIN.E( I ))2 

PXY. E(I)  - 2.*SANG*CANG*(1 . /AHAJ.E( I )2  - 1/AMIN. E(I )?) 

The  following  SIMSCRJPT  1 1. 5 code  segment  then  computes  the  tree  height 

at  an  arbitrary  point  X,Y  assuming  that  there  are  NCvtLS  cover  ellipses  with 

' \ 

global  computing  parameters  as  defined  Above. 


TREE 

.HT 

■ 0. 

I - 

1 TO  NOVELS  DO 

LET 

XS 

■ X-XC.ECI ) 

LET 

YS 

= Y-YC.E(I) 

LET 

QI 

« PXX.E(I)*XS**2  ♦ PYY.E( I )*YS**2 

♦ PXY.E(I)*XS*YS 

IF 

QI 

GE  1 . CYCLE 

ELSE 

IF 

HT.E(I)  GT  TREE.HT  LET  TREE.HT 

ALWAYS 

LOOP 

0.  Efficiency 

As  Is  the  case  for  terrain  hills,  greater  efficiency  can  be  obtained 


- 
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by  Uniting  the  number  of  cover  ellipses  checked  at  each  X,Y  to  those 
found  In  a grid  square.  A forest  preprocessor  program  called  TREE. LIST 
operates  In  a fashion  exactly  analogous  to  HILL. LIST  to  develop  an  array 
LIST.C(iX,lY,L)  of  cover  ellipse  numbers  In  grid  square  IX, IY.  Again  this 
Is  stored  as  a SIKSCRIPT  ragged  array.  The  resulting  more  efficient  code 
segment  differs  from  the  above  code  just  as  In  the  ELEV  program.  Since  It 
appears  In  the  LOS  program  we  will  not  detail  It  further  here  beyond  de- 
fining the  LIST.C  array  contents.  For  given  IX, IY 

LIST.C(IX,IY,1)  ■ number  of  forest  ellipses  overlapping 
grid  square  IX.IY. 
and 

LIST.C(IX,IY,2),...etc.  give  the  ellipse  numbers  If 
LIST. C( IX,IY,1 ) > 1 . 

Documentation  of  the  TREE. LIST  program  appears  In  Section  VII. 

V.  Overview  of  LOS  Modelling  In  STAR 
An  Important  computation  for  any  high  resolution  combat  simulation  Is 
the  line  of  sight  (LOS)  routine.  The  basic  problem  Is  the  following:  Given 
an  observer  (A)  and  a potential  target  (B),  what  part  (If  any)  of  the  target 
can  be  seen  by  the  observer?  LOS  Is  a purely  geometric  computation  in  that 
we  assume  perfect  visibility.  Degraded  visibility  conditions  are  Incorporated 
In  other  modules  of  the  STAR  model.  The  result  of  the  LOS  computation  Is  a 
percent  of  the  vertical  height  of  the  target  B visible  to  the  observer  A. 
This  percent  visible  Is  used  In  various  ways  In  other  parts  of  the  STAR  model. 

The  LOS  computations  are  rather  complex,  but  conceptually  they  are 
based  on  a simple  procedure:  "Find  the  lowest  sight  line  from  A over  the 
terrain  (and  forests).  Extend  this  line  to  B's  location,  and  compare  its 
extrapolated  height  to  B's  elevation.  Thus  compute  % visible."  This 
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procedure  Is  Illustrated  In  figure  13  for  two  cases.  In  the  first,  a terrain 
hill  blocks  LOS  completely  so  0%  of  the  target  is  visible.  In  the  second 
case  the  lowest  sight  line  intersects  the  target  so  that  40&  of  its  vertical 
height  can  be  seen  by  A while  the  other  60X  Is  covered. 

This  simple  concept  is  rather  difficult  to  compute  because 

1.  any  one  of  the  hills  or  forests  which  lie  between  A and 
B may  be  the  feature  which  determines  the  lowest  sight 
line.  Thus  we  have  to  repeat  computations  for  each  such 
feature,  and 

2.  the  equation  for  the  lowest  sight  line  over  a hill  does  not 
have  a closed  form  solution. 

The  following  observations  tend  to  lessen  the  above  computational  difficulties: 

3.  Since  A and  B are  given  locations,  only  hills  and 
forests  which  intersect  the  line  joining  A and  B need  to 
be  considered.  Features  which  are  remote  from  this  line 
can  be  Ignored  totally. 

4.  If  LOS  does  not  exist  at  all,  (01  visible)  then  as  soon  as 

a hill  or  forest  has  been  found  which  proves  that  the  target 
is  01  visible,  the  computation  can  stop.  All  other  hills 
and  forests  can  be  ignored. 

Point  4.  is  particularly  significant  because  on  a typical  battlefield 
a large  fraction  of  the  potential  observer-target  pairs  will  not  have  line 
of  sight.  Since  the  lowest  sight  line  computation  is  difficult,  the  STAR  LOS 
routine  first  performs  a number  of  simpler  tests  trying  to  prove  either  that 
% visible  a 0 , or  that  a particular  hill  Is  so  small  that  it  could  not 
possibly  Interrupt  LOS  between  A and  B.  In  either  case  the  difficult 
lowest  sight  line  computation  is  avoided,  and  in  the  first  case  all  further 
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computations  are  unnecessary.  Only  in  cases  where  % visible  is  non  zero 
do  we  have  to  do  the  difficult  computations  for  many  hills. 

In  this  section  we  briefly  discuss  the  various  tests  in  the  LOS  routine. 
The  discussion  will  be  in  terms  of  the  geometry  of  the  LOS  situation.  The 
derivation  of  the  associated  mathematical  computations  is  deferred  to 
Section  VI  of  this  report.  LOS  assumes  the  following  data  is  available  for 
each  of  A and  B: 

XA,YA(XB,YB)  The  X and  Y coordinates  on  the 
battlefield  for  each  of  A and  B . 

TMACA, (TMACB)  The  macro  terrain  elevation  computed 
from  the  ELEV  routine. 

TMICA  (TMICB)  A micro  terrain  offset  + or  - from 

the  macro  terrain.  This  offset  is  used  in 
several  ways: 

i)  For  stationary  defenders  in  defilade, 
negative  micro  terrain  offsets  simulate 
the  defilade  position. 

ii)  For  elements  in  the  open,  a random  micro 
terrain  choice  can  be  used  to  simulate 
terrain  irregularities  that  are  too  small 
to  explicitly  model  in  the  ELEV  routine, 

iii)  Large  positive  micro  terrain  elevations 
are  used  to  simulate  aircraft  flying 
above  the  terrain. 

SIZEA  (SIZEB)  The  vertical  height  dimension  of  each  element. 

The  percent  visible  is  taken  as  a fraction  of 
this  size. 
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LAGA  (LAGB) 


Indicator  variables  for  air/ground 
0 =?>  the  element  is  on  the  ground 
1=^  the  element  is  an  aircraft 


LATOB  (LBTOA)  Indicator  variables  for  one  or  two 
way  LOS  calls. 

LATOB  = 1 compute  LOS  from  A to  B 
(yielding  VISFRB) 

0— >do  not  compute  A to  B 

Since  many  of  the  situations  in  which  LOS  is  computed  are  two-way  situations, 
(Eg.  A is  trying  to  detect  B,  but  B is  also  trying  to  detect  A)  the 
routine  includes  the  option  to  do  all  calculations  in  both  directions,  thus 
saving  repetition  of  many  common  computations.  The  result  of  the  LOS  compu- 
tations Is  provided  in  the  two  variables 

VISFRA  (VISFRB)  The  fraction  of  SIZEA  (SIZEB)  which 
can  be  seen  by  B(A) 

To  compute  LOS  we  evaluate  terrain  and  forest  height  at  a number  of 
points  along  the  line  segment  between  A and  B.  The  various  tests  are 
sequenced  in  the  same  order  in  the  computer  code  as  in  the  discussion  which 
follows.  For  this  discussion  we  will  consider  a one-way  LOS  CALL  with  A as 
observer  and  B as  target  (LAT0B=1,  LBTOA  = 0).  The  two  way  computation  is 
exactly  analogous. 

A.  Initialize 

To  begin  with,  the  LOS  routine  sets  the  visible  fraction  VISFRB=1.0  . 
As  we  proceed  through  the  routine,  various  tests  may  decrease  the  visible 
fraction,  but  none  will  ever  increase  it.  If  VISFRB  ever  decreases  to  0, 
then  we  immediately  terminate  the  computations  and  return  from  the  routine. 


* 

f 
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Throughout  the  computations  we  assume  that  the  bottom  of  A Is  at 
TMACA  ♦ TMICA  and  that  the  top  of  A is  at  ZA  = TMACA+TMICA+SIZEA  and 


similarly  for  B (see  Figure  14).  The  observation  device  is  assumed  to  be  at 
the  top  of  A.  If  the  micro  terrain  value  TMICB  Is  negative,  then  it  is 
assumed  that  the  bottom  part  of  the  target  Is  not  visible,  so  in  this  case 
VISFRB  Is  Immediately  decreased  to 

VISFRB  - 1.  ♦ (19) 

which  Is  less  than  1 when  TMICB  < 0 . 

B.  Grid  Square  List 

The  line  segment  between  A and  B's  positions  crosses  over  some 
of  the  terrain  grid  squares.  Thus  the  only  hills  and  forests  we  need  to  con- 
sider are  those  associated  with  these  grid  squares.  The  routine  develops  a 
list  of  the  NGRSQ  grid  squares  crossed  in  the  arrays  IGX(K)  and  IGY(K)  for 
K=1 .... .NGRSQ. 

C.  Forest  Ellipse  List 

The  line  segment  between  A and  B may  intersect  some  of  the 
forest  cover  ellipses.  If  so,  we  must  know  where  the  intersection  points  lie. 
The  routine  develops  a list  of  these  ellipse  intersections  by  considering, 
for  each  grid  square  only  the  ellipses  (in  LIST.C)  which  appear  in  the  grid, 
and  then  testing  each  for  intersection  with  the  A to  B line.  The  inter- 
section points  (if  they  exist)  are  called  SI  and  S2  and  are  stored  along 
with  the  ellipse  numbers  in  arrays  IEL(K),  CS1(K)  and  CS2(K)  for  K=1 .... ,NELS. 
Although  an  ellipse  may  appear  in  several  adjacent  grid  squares,  it  will  be 
stored  on  the  list  at  most  one  time. 

D.  LOS  Test  at  Forest  Boundaries 

One  point  at  which  lack  of  LOS  can  often  be  proved  is  a forest 


j 
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boundary,  so  as  the  above  list  of  forest  ellipses  is  being  developed  LOS  is 
tested  In  the  following  ways: 


1.  If  A or  B Is  a ground  element  and  Is  inside  a 

forest  ellipse,  VISFRB=0.  (no  visibility  through 
forests) 

2.  If  A or  B is  an  air  element  but  the  corresponding 

TMIC  value  is  less  than  forest  height  (e.g.  a 

hell  coper  hovering  among  the  trees,  but  below  the 
tree  line)  then  VISFRB=0  . 

3.  If  SI  (or  S2)  is  between  A and  B then  the  macro 
terrain  elevation  is  added  to  the  tree  height  at  SI 
(or  S2)  and  the  combined  elevation  is  tested  to  see 
if  it  interrupts  or  decreases  the  percent  visible. 

This  situation  is  diagrammed  in  cross-section  in 
Figure  15.  Subroutine  TREE. CHECK  performs  these 
computations. 

E.  HILLTOPS 

Next  we  consider  all  macro  terrain  hills  that  might  lie  between 
A and  B.  The  list  of  such  hills  is  developed  by  scanning  LIST.H  for  each 
grid  square  crossed  by  the  A to  B line.  Each  such  hill  is  processed  one 
at  a time  from  this  point  through  to  the  end  of  tie  computations.  Although  a 
hill  may  appear  in  several  grid  squares,  each  hill  is  considered  only  once. 

We  first  compute  a point  W on  the  A to  B line  which  is  the  top  of  the 
hill  cross-section  cut  by  that  line.  The  following  tests  are  then  done  on  W. 

1.  If  W is  far  from  being  between  A and  B,  then 
this  hill  is  irrelevant;  go  on  to  the  next  hill 
(see  Figure  16,  part  1). 
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2.  If  the  hill  height  at  W Is  negligible  (1,*  the 
cross  section  Is  far  from  the  hill  center)  then 
this  hill  Is  Irrelevant;  go  on  to  the  next  hill. 

(See  Figure  16,  part  2). 

3.  Test  the  hill  height  at  W plus  the  tree  height 

at  W (If  any)  to  see  If  LOS  Is  totally  Interrupted. 

(See  Figure  16,  part  3). 

4.  If  hill  height  at  W plus  thr  tree  height  Is  below 
the  bottoms  of  both  A and  B,  then  this  hilltop 
cannot  affect  LOS,  go  on  to  the  next  hill. 

(See  Figure  16,  part  4). 

F.  Lowest  Sight  Line 

If  no  previous  computations  have  reduced  VISFRB  to  0.,  then  LOS 
probably  exists,  and  we  must  compute  the  lowest  sight  line  over  the  hill  to 
compute  the  actual  percent  visible.  It  is  not  enough  to  consider  the  hilltop, 
since  for  cases  where  A Is  significantly  above  or  below  the  hilltop,  the 
lowest  sight  line  will  graze  the  hilltop  away  from  W . (See  Figure  17).  We 
denote  by  V the  point  at  which  the  lowest  sight  line  is  tangent  to  the  hill. 
Once  V Is  found,  then  the  lowest  sight  line  is  extrapolated  to  B's  position 
for  the  computation  of  VISFRB  (as  In  Figure  13).  Any  forest  coverage  at  V Is, 
of  course,  also  considered. 

When  these  computations  have  been  done  for  each  hill,  with  VISFRB 
possibly  decreased  at  each  step,  the  final  resulting  VISFRB  value  Is  returned 
as  the  result  of  the  LOS  routine. 

VI . Details  of  the  Line  of  Sight  Computations 
In  this  section  we  derive  the  mathematical  basis  for  the  computations 
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In  the  LOS  routine  and  examine  the  SIMSCRIPT  LOS  code  line  by  line.  This 
section  is  the  most  technical  of  the  report,  and  the  following  sections  can  be 
read  without  comprehending  all  of  this  material.  Thorough  understanding  of 
this  section  is,  however,  required  for  anyone  who  wants  to  adjust  or  modify 


the  STAR  LOS  procedure.  The  computations  will  be  detailed  in  the  same  order 
as  they  were  discussed  In  Section  V,  and  in  essentially  the  same  order  as  they 
appear  In  the  code.  In  the  process,  the  complete  LOS  code  will  be  listed  and 
briefly  annotated.  Throughout  the  discussion  we  use  the  same  notation  for 
variables  and  constants  as  in  the  LOS  program  with  the  following  exception: 

To  speed  execution  of  the  frequently  called  LOS  program,  its  many  local 
variables  have  all  been  made  global.  To  assure  unique  names  for  these  global 
variables,  and  to  aid  in  identifying  them,  each  such  variable  name  has  been 
suffixed  with  .LS,  so  that  TMACA  in  our  discussions  will  appear  as 
TMACA.LS  In  the  computer  code. 

The  mathematics  of  the  LOS  routine  depends  crucially  on  the  following 

observation.  Consider  the  elevation  of  a single  hill  I given  by 

Z(X,Y)  = PEAK.H(I)  - HT.H(I)  + HT.H(I)  * exp  [QT(X,Y)] 

1 (20) 

= Cl  + C2  * exp  [Qr(X,Y)] 

where  Qj  Is  a negative  definite  quadratic  function  of  X and  Y in  2-space. 
Parameterize  the  straight  line  In  X,Y  space  between  an  observer  A and  a 
target  B as  follows: 

X(S)  = XA  + S*(XB  - XA) 

(21) 

Y(S)  * YA  + S*(YB  - YA) 

so  that  when  S=0  we  are  at  A's  location,  and  when  S=1  we  are  at  B's 
location.  Then  consider  the  elevation  of  hill  I as  a function  of  S along 
the  A to  B line.  The  observation  is  that  the  elevation  is  given  by  a 
function  of  the  same  form  as  (20)  except  now  with  a single  variable  S . 


Z(S)  = Cl  + C2  exp  [P(S)]  (22) 

where  P(S)  Is  a negative  definite  quadratic  function  of  the  parameter  S . 

To  obtain  the  coefficients  o * P(S)  we  plug  X(S)  and  Y(S)  from 
(21)  into  (20).  For  notational  simplicity  let  XC,  YC,  PXX,  PYY,  PXY  denote 
the  computational  parameters  for  hill  I,  then  we  have 
Qj(X,Y)  = Qj(X(S),Y(S))  = 

• PXX*(X(S)-XC)2  + PYY*(Y(S)-YC)2  + PXY*(X(S)-XC)*(Y(S)-YC) 

= PXX*(XA+S*XBA-XC)2  + PYY* ( YA+S*YBA- YC ) 2 
+ PXY*(XA+S*XBA-XC)*( YA+S*YBA- YC ) 

(where  XBA  = XB  - XA  and  YBA  = YB  - YA) 

= PXX*(RX+S*XBA)2  + PYY*(RY+S*YBA)2 

+ PXY*(RX+S*XBA)*(RY+S*YBA)  ! 

(where  RX  = XA  - XC  and  RY  = YA  - YC) 

■ S2*[PXX  * XBA2  + PYY  * YBA2  + PXY  * XBA  * YBA] 

+ S*[2PXX  * XBA  * RX  + 2PYY  * YBA  * RY 

+ PXY  * (XBA  * RY  + YBA  * RX)] 

+[PXX  * RX2  + PYY  * RY2  + PXY  * RX  * RY] 

= S2  * GQ  + S * FQ  + EQ  defining  the  quadratic  coefficients 
GQ.FQ.EQ  and  giving  the  formula  for  P(S).  (23) 

Given  this  equation,  all  future  LOS  computations  can  deal  with  the 
single  dimensional  computation  along  the  A to  B line  parameterized  by  S . 

We  begin  the  discussion  with  a simple  subroutine  which  is  called  from 
LOS  in  several  places.  Routine  KOVER  abstracts  the  LOS  percent  visible  com- 
putation to  its  bare  essence.  Consider  parameterizing  an  observer-to-target 
line  in  XY  space  with  a distance  parameter  S so  that  S=0  represents  the 
observer's  position  and  S=1  represents  the  target's  location.  Assume: 

(See  Figure  18) 
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1.  The  observation  device  has  elevation  Z0. 

2.  The  target  has  a macro-terrain  elevation  TMACT,  a micro  terrain 
offset  TMICT  and  a size  of  SIZET.  Thus  its  top  is  at 

ZT  •=  TMACT  + TMICT  + SIZET. 

3.  At  some  point  S between  observer  and  target  (0  £ S £ 1)  there  is 
an  obstruction  to  line  of  sight  whose  top  has  elevation  HTS. 

4.  Other  previous  computations  have  established  a visible  fraction 
of  VISFIN  due  to  other  obstructions. 

To  determine:  Does  the  obstruction  at  S reduce  the  percent  visible? 

Let  ZS  - Z0  + S*(ZT  - ZS)  (24) 

If  HTS  >_  ZS,  then  the  obstacle  totally  blanks  LOS  and  the  resulting  percent 
visible  Is  VISF0UT=0. 

If  HTS  < ZS, 

Let  EVIST  = max  (TMACT,  Zf>  + (HTS  - Z0)/S)  (25) 

giving  the  elevation  of  the  lowest  point  on  the  target  which  can  be  seen, 
if  EVIST  >_  ZT  then  LOS  is  blocked,  and  otherwise 

VISFOUT  = min  (VISFIN,  ( ZT-E VI ST ) /S I ZET ) (26) 

« 

possibly  decreasing  the  previously  computed  visible  fraction  VISFIN.  Given 

the  above  formulas,  the  code  for  routine  KOVER  should  be  self-explanatory.  All 
variables  in  the  routine  are  local.  KOVER  is  only  called  from  LOS  or  from 
other  routines  which  LOS  calls.  For  the  KOVER  program  see  Listing  3. 

We  now  consider  the  various  segments  of  the  LOS  routine  in  Listing  4. 

A.  Initialize 

In  addition  to  initializing  VISFRA  and  VISFRB,  this  segment  of 
the  code  computes  some  temporary  values  which  are  used  repeatedly  later  in  the 
code.  KTREP  is  a global  counter  which  is  initialized  to  -INF.C  in  RES. TERR 
and  increases  by  one  on  each  call  to  LOS.  It  is  used  to  avoid  processing  a 
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given  hill  or  forest  ellipse  more  than  one  time.  Note  in  lines  9 and  10 
that  if  the  micro  terrain  offset  is  negative,  then  part  of  the  target  is 
assumed  masked  by  the  terrain,  so  the  visible  fraction  starts  out  at 

less  than  1 .0  . 

B.  Grid  Square  List 

The  grid  square  list  is  developed  starting  at  grid  square  IX, IY 
which  contains  point  B (lines  29-34).  This  grid  square  goes  on  the  list  as 

IGX(l),  I GY ( 1 ) . ISGX,  ISGY  are  set  to  + 1 depending  on  whether  the  grid 

square  indices  IX,  IY  increase  or  decrease  as  we  move  from  B to  A. 

(Lines  20-27).  XSTEP  and  YSTEP  are  initialized  (in  lines  35-36)  to  the 
fraction  of  the  distances  XBA  and  YBA  we  have  moved  as  we  map  out  the  grid 

squares  from  B to  A . XINC  and  YINC  are  set  to  the  increment  in  XSTEP 

and  YSTEP  which  is  added  when  we  move  across  one  grid  square  in  the  X or 
Y directions  (Lines  20-27). 

The  loop  from  lines  37  to  47  sets  XSTEP  and  YSTEP  to  the  next 
X-grid  intersection  and  the  next  Y-grid  intersection.  The  smaller  of  XSTEP 
and  YSTEP  signals  which  intersection  we  hit  first  and  hence  which  of  IX,  IY 
gets  incremented  by  ISGX  or  ISGY.  After  each  increment  a new  grid  square 
is  recorded  in  IGX,  IGY.  When  both  XSTEP  and  YSTEP  are  greater  than  1.0 
we  have  passed  A,  and  the  loop  terminates  with  NGRSQ  grid  square  indices 
recorded  in  the  arrays  IGX,  IGY. 

C.  Forest  Ellipse  List  and  D.  LOS  Test  at  Forest  Boundaries 

If  there  are  no  forest  ellipses  on  the  battlefield  (NCVELS=0) 

then  this  segment  is  bypassed.  (Line  52).  Otherwise  we  loop  over  all  grid 
squares  crossed,  and  for  each  grid  square  loop  over  all  forest  ellipses,  IC, 
which  intersect  the  grid  square  (lines  53-58).  Lines  59-60  ensure  that  each 
ellipse  is  accessed  at  most  one  time. 
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For  each  such  ellipse,  we  want  to  know  whether  or  not  the  A to  B 


line  intersects  the  ellipse,  and  if  so,  where.  The  ellipse  boundary  equation 
is  stored  as 

PXX*(X-XC)2  + PYY*(Y-YC)2  + PXY*( X-XC)* (Y-YC)  = 1 (27) 

and  X,Y  on  the  A to  B line  are  given  parametrically  as  in  (21). 

Solving  (21)  and  (27)  simultaneously  for  S gives  the  intersection  points 
(if  any)  as  follows: 

PXX*(XA-XC  + $*XBA)2  * PYY*YA-YC  ♦ X*YBA)2  ♦ 

PXY*(XA-XC  + $*SBA)*CYA-YC  + S*YBA)  - 1 

o ? (28) 

PXX*(RX  + S*XBAT  + PYY*(RY  + S*YBA;  + 

PXY*(RX  + S*XBA)*(RY  + S*YBA)  * 1 

where  RX  = XA  - XC  and  RY  = YA  - YC  . 

This  simplifies  as  in  (23)  to  the  simple  quadratic  equation  in  S; 

AA  = PXX*XBA2  + PYY*YBA2  + PXY*XBA*YBA 

BB  = 2PXX*XBA*RX  + 2PYY*YBA*RY  + PXY*(XBA*RY  + YBA*RX)  (30) 

CC  = PXX*RX2  + PYY*RY2  + PXY*RX*RY  - 1 .0 
Applying  the  quadratic  formula  gives  the  intersection  points. 

Cl  _ -BB-  V BB2-4.*AA*CC 

S1 — - <31) 


S2 


-BB+ 


-4.*AA*CC 


(32) 


(with  S,1  ^ S2)  If  the  radical  is  positive,  and  no  intersection  if  the  radical 
is  negative.  Lines  61-71  perform  this  computation. 

If  SI  and  S2  exist  then  we  need  to  consider  whether  the  tree  height 
at  the  forest  boundaries  SI,  S 2 Interrupts  LOS.  Oepending  on  whether  one 
or  both  of  SI  and  S2  He  between  A and  B (that  is  between  S=0  and 
S=l),  and  depending  on  whether  A and  B are  air  or  ground  platforms,  it 

may  be  necessary  to  compute  macro  terrain  elevation  plus  tree  height  at 
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either  or  both  of  SI  and  S2  and  see  if  the  resulting  obstacle  decreases 
the  percent  visible.  Subroutine  TREE. CHECK  does  the  terrain  computation  and 
LOS  check  (using  KOVER),  and  will  be  documented  following  the  LOS  code. 

Lines  (72-94)  enumerate  the  possible  cases  which  determine  if  TREE. CHECK  is 
to  be  called  at  SI  and/or  S2  . 

Finally,  for  this  segment,  the  ellipse  number  IC,  and  the  SI,  S2 
values  for  any  ellipse  which  intersects  the  A-to-B  line  are  saved  in  arrays 
IEL,  CS1,  CS2  for  later  use  in  the  LOS  procedure.  (Lines  95-98). 

E.  Hilltops 

Next  the  macro-terrain  hills  lying  between  A and  B must  be 
checked  to  see  if  they  interrupt  LOS.  As  indicated  in  Section  V,  we  loop 
over  all  NGRSQ  grid  squares  and  for  each,  consult  the  LIST.H  array  to  qet  the 
hill  numbers  I to  be  considered.  (Lines  105-111).  As  for  the  forest 
ellipses,  if  a hill  has  already  been  considered  in  another  grid  square,  then 
it  is  skipped  over  (Lines  112-113). 

The  only  interesting  computation  in  this  section  is  finding  the 
location  W and  height  HHW  of  the  hilltop.  In  general  W will  not  be  at  the 
hill  center,  since  the  A to  B line  need  not  pass  through  the  center.  Thus 
HHW  is  usually  somewhat  less  than  PEAK.H(I).  The  hilltop  location  is  easy  to 
compute  because  of  equations  (22)  and  (23)  which  give  the  hill's  elevation 
along  the  A to  B line  as 

Z(S)  = PEAK . H ( I ) -HT . H ( I ) + HT.H( I )*exp[GQ*S2  + FQ*S  + EQ]  (33) 
The  hilltop  must  occur  where  dZ/dS  = 0,  so  differentiating  (33)  and  setting 
the  result  to  zero  gives 

= HT.H(I)*exp[GQ*S2+FQ*S+EQ]*[2GQ*S+FQ]  =0  (34) 

if  and  only  if 

[2GQ*S  + FQ]  = 0 (35) 

- 55  - 


or 


W = S = - FQ/(2*GQ)  (36) 

Lines  114-124  compute  GQ,  FQ,  EQ  as  In  (23)  and  W as  in  (36).  The  hill 
height  HHW  at  W is  obtained  by  substituting  W back  into  (33)  giving 

HHW=Z(W)=PEAK.H( I ) - HT.H(I)  + HT.H(I)*exp[EQ-FQ2/(4*GQ))]  (37) 

in  (lines  125-128).  The  tests  of  Figure  16  are  performed  in  lines  129-141, 
and  in  particular,  lines  133-139  loop  through  the  previously  saved  forest 
ellipses  to  compute  the  tree  height  (if  any)  at  W . 

F.  Lowest  Sight  Line 

The  lowest  sight  line  procedure  to  compute  the  tangency  point  V , 
is  an  iterative  procedure  using  the  Newton-Raphson  equation  solving  pro- 
cedure. This  is  necessary  because  the  equations  describing  the  lowest  sight 
line  do  not  have  a closed  form  solution.  In  the  LOS  code  this  computation  is 
delegated  to  subroutine  NEWTON  (to  be  described  shortly).  The  NEWTON  procedure 
uses  slightly  different  parameters  depending  on  the  direction  of  the  computa- 
tion A to  B or  B to  A . Both  calls  are  set  up  in  lines  144-153  of 
the  LOS  code.  The  details  of  the  NEWTON  iteration  will  be  documented  along 
with  routine  NEWTON. 

Finally,  the  LOS  routine  returns  to  the  calling  program  either  with  the 
computed  visible  fractions  (on  line  157)  after  all  hills  have  been  checked  in 
all  NGRSQ  grid  squares,  or  with  VISFRA  = VISFRB  = 0.0  if  no  LOS  exists  (line 
159).  This  completes  the  detailed  documentation  of  the  LOS  routine. 

G.  Routine  TREE. CHECK 

Given  element  A at  S=0,  element  B at  S=l,  and  a forest 


boundary  at  a point  SS=S1  or  S2  with  0 <_  SS  <_  1 , the  TREE. CHECK  routine 
tests  whether  the  forest  edge  Interrupts  LOS  from  A to  B or  from  B to  A. 
The  routine  has  no  local  variables,  sharing  the  LS  global  variables  with  LOS. 
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The  TREE. CHECK  Code  is  given  in  Listing  5 
Line  2 computes  the  X and  Y coordinates  of  the  forest 
boundary  at  SS  using  equation  (21). 

Line  3 calls  ELEV  to  get  the  macro  terrain  elevation  HTS  at  SS 
Line  4 adds  to  HTS  the  tree  height  for  this  forest  ellipse 

Lines  6-12  set  up  the  A to  B and  B to  A calls  to  routine 
KOVER  to  see  if  percent  visible  is  decreased  by  the 
forest  obstacle. 

H.  Routine  NEWTON 

Given  the  equation  (33)  for  a hill  parameterized  in  S,  and  given 
an  observer^  location  at  S=0  with  observation  device  at  elevation  ZA,  we 
wish  to  compute  the  point  S=V  at  which  the  lowest  sight  line  is  tangent  to 
the  hill  (if  such  a point  exists),  (See  Figure  17).  The  elevation  of  the 
sight  line  at  any  S is  given  by  the  equation 

L(S)  = ZA  + S*SL0PE  (38) 

While  the  hill  elevation  is,  as  before, 

Z(S)  = PEAK.H-HT.H  + HT.H*exp  [GQ*S2  + FQ  * S + EQ]  (39) 

pTs) 

Tangency  at  V requires  that  the  elevations  be  equal, 

L ( V)  = Z (V)  (40) 

and  also  that  the  slopes  be  equal  at  V . 

L ' (V)  = Z 1 ( V ) (41) 

Unfortunately  this  set  of  equations  does  not  have  any  closed 
form  solution  known  to  this  author.  Instead  we  apply  the  Newton  Raphson 
iteration  to  solve  for  V in  the  following  fashion: 
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Z(V)  = L(V)  = ZA  + L'(V)*V 

by  (40)  and  since  L ( S ) is  linear.  Thus,  using  (41), 

Z(V)  = ZA  + Z'(V)*V 

Also,  from  (34) 

Z • (V)  = HT.H*exp[P(V)]*P’(V) 

Plugging  (39)  and  (44)  into  (43)  gives  an  equation  for  V: 

PEAK.H  - HT.H  + HT.H*exp[GQ*V2  + FQ*V  + EQ]  = 

ZA  + V*HT.H*exp[GQ*V2  + FQ*V  + EQ]*[2GQ*V+FQ] 


(42) 

(43) 

(44) 


(45) 


Or,  collecting  terms 


FCN (V)£  HT.H*exp[GQ*V2+FQ*V+EQ]*(2GQ*V2+F*V-l ) 


(46 


+ ZA  + HT.H  - PEAK.H  = 0 
To  find  the  V which  makes  FCN(V)  = 0,  we  use  the  iterative  Newton  Raphson 
procedure,  starting  with  the  reasonable  guess  = W and  repeatedly  applying 
the  iteration  formula 


Vi  = vk 


FCN(Vk) 

DFCN(Vk) 


(47) 


where  DFCN  is  the  derivative  of  FCN.  If  the  iteration  leads  to  a V value 
which  is  remote  from  the  location  of  A at  S=0  and  B at  S=l,  then  the 
iteration  is  abandoned,  and  we  cycle  to  the  next  hill.  If  the  iteration 
converges  to  a V fo  which 

Z(V) « ZA  + Z'(V)*V  (48) 

then  this  V is  tested  for  its  effect  on  LOS  by  calling  the  KOVER  routine. 

The  actual  code  which  follows  in  Listing  6 is  slightly  more  involved 
than  this  description  since  it  must  be  able  to  compute  LOS  in  both  directions. 
This  is  accomplished  using  the  parameter  VSUB  which  is  0 for  an  A to  B 
call  and  1 for  a B to  A call. 
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I . Interface  of  LOS  with  STAR 

To  use  the  LOS  routine  the  STAR  model  must  set  14  parameters: 

(all  global  variables  with  the  .LS  suffix) 

XA,  YA,  TMACA,  TMICA,  SIZEA,  LAGA,  LATOB 

XB,  YB,  TMACB,  TMICB,  SIZEB,  LAGB,  LBTOA 

and  CALL  LOS.  On  return  VISFRA  and  YISFRB  have  the  fraction  of  A visible 

to  B and  the  fraction  of  B visible  to  A respectively.  If  a one-way 

call  was  requested  (say  LBTOA  = 0)  then  one  of  the  visible  fractions  (VISFRA 

in  this  case)  will  contain  a meaningless  value  and  should  be  ignored. 

Since  the  most  frequent  use  for  LOS  is  between  entities  in  the  sim- 

(see  Listing  7) . 

ulation,  a driver  routine  SIGHT  has  been  written./  Given  the  pointers  A and 
B to  the  two  entities,  SIGHT  sets  up  the  14  global  variables  from  the 
appropriate  entity  attributes,  and  calls  LOS.  Note  that  this  routine  assumes 
both  A and  B are  ground  elements.  The  generalization  to  air  elements  is 

I ' 

tr  vial.  On  return,  the  global  variable  PCT.VIS  is  set  to  VISFRB  scaled  by 
the  size  of  the  target  B relative  to  the  size  of  the  M60  tank  because  the  de- 
tection model  in  STAR  is  based  on  the  DYNTACS  field  experiments  which  used 
the  M60  as  a target. 

VII . Terrain  Preprocessors  and  Data  Input 
The  efficient  computation  of  the  STAR  macro-terrain  and  forest 
representations  requires  prior  execution  of  the  HILL. LIST  and  TREE. LIST  pro- 
grams and  the  routine  RES. TERR  which  reads  terrain  and  forest  data  into  the 
simulation. 

A.  HILL. LIST 

As  indicated  in  Section  III,  the  HILL. LIST  program  prepro- 
cesses the  macro  terrain  hill  data  to  improve  execution  efficiency.  The 
preprocessing  consists  of  the  following  steps. 
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1.  Define  the  battlefield  and  an  array  of  grid 
squares  covering  It. 

2.  Scan  each  grid  square  to  build  a list  of  which 
hills  actually  contribute  to  the  macro  terrain 
In  that  square. 

3.  For  each  grid  square,  sort  the  list  of  hills 

so  that  the  biggest  are  first  on  the  list.  Thus 
If  a hill  interrupts  LOS  we  are  likely  to  find 
It  sooner. 

4.  For  each  hill,  find  the  lowest  point  on  that  hill 
that  actually  contributes  to  the  macro  terrain,  and 
thus  (possibly)  improve  the  cutoff  value  for  the  hill. 

5.  Output  a deck  of  data  cards  with  hill  parameters 
and  the  LIST.H  array  ready  to  be  read  into  STAR  by 
routine  RES. TERR. 

The  Code  is  given  in  Listing  8. 


Input  to  HILL. LIST  is  SIMSCRIPT  free  format  from  cards.  The  following  values 


are  required  (in  order):  (Lines  5-6) 


NGRIDX 

NGRIDY 


4 

number  of  grid  squares  in 
' X and  Y directions 

, 


GSIZE 

SPACING 

X. LO.BDRY 

Y. LO.BDRY 
NHILLS 


size  in  meters  of  each  grid 
square  (eg.  1000.) 

increment  in  meters  for  scanning 
each  grid  square  (eg.  50.) 

battlefield  coordinates  in  meters 
of  the  southwest  corner  origin 
of  the  grid  square  system 

number  of  hills  to  be  used  on  the 
battlefield. 


n 


‘ 


1 


* ■ 
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Next  the  fitting  parameters  for  each  hill  are  input:  (line  17-18) 

XC.H,  YC.H,  PEAK.H,  ANG.H,  ECC.H,  SPRD.H,  HT.H. 

The  current  program  does  not  read  CUT.H  since  cutoffs  were  not  used  in  the 
original  battlefield  coding.  These  could  easily  be  added.  The  fitting 
parameters  are  converted  to  computing  parameters  in  lines  19-27. 

Finally,  BASE  values  for  each  grid  square  are  entered.  The  current 
code  has  these  hard  wired  since  only  3 distinct  base  values  were  used  on 
the  original  battlefield.  A card  read  could  easily  replace  this.  (Lines 
30-35). 

Lines  37-59  perform  the  battlefield  scan,  computing  elevation  every 
SPACING  meters  and  recording  which  hill  was  responsible  for  the  final  ele- 
vation. This  hill  number  is  entered  in  a temporary  array  TEMP  if  it  is  not 
already  there.  (Lines  48-55).  Array  ZHILL  records  the  contribution  of  each 
hill  to  this  grid  square  for  the  upcoming  sort.  h 

Lines  60-78  sort  the  TEMP  array  and  store  it  in  the  appropriate  column 
of  LIST.H.  The  rest  of  the  program  prints  and  punches  the  output  data  deck. 

This  deck  is  ready  for  immediate  input  to  the  simulation  via  routine  RES .TERR  . 

B.  TREE .LIST 

in  Listing  9 

The  TREE. LIST  program/\has  exactly  the  same  structure  as  the  HILL. LIST 
program.  It  is  presented  here  with  essentially  no  further  comment.  The  data 
output  deck  from  TREE. LIST  is  ready  to  be  read  by  RES  .TERR  except  that  the 
first  card  (battlefield  definition)  which  duplicates  the  first  card  from 
HILL. LIST  should  be  discarded. 

C.  RES. TERR  Routine 

RES. TERR  is  a routine  called  from  the  STAR  MAIN  program  which 
initializes  the  following  data: 
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1.  Battlefield  definition  and  grid  square  system 

2.  Macro-terrain  hill  computing  parameters 

3.  LIST.H  array 

4.  Forest  ellipse  computing  parameters 

5.  LIST.C  array 

Input  to  RES. TERR  consists  of 

1.  The  punched  output  deck  from  HILL. LIST  followed  by 

2.  Either  a single  card  with  a 0 (zero)  if  there  are 
no  forest  ellipses  or  else  the  punched  output  deck 
from  TREE. LIST  minus  its  first  card  (which  duplicates 
the  battlefield  definition  card  from  HILL. LIST) 

RES. TERR  is  set  up  to  dynamically  reserve  and  dimension  the  various  arrays 
so  that  no  more  core  is  required  than  is  absolutely  necessary.  For  example, 

if  there  are  no  forest  ellipses,  then  no  data  arrays  are  reserved  for 

ellipse  parameters  or  for  LIST.C.  This  first  implementation  of  RES. TERR 
makes  no  attempt  to  pack  data  values,  although  if  required  some  packing  should 
be  possible.  The  program  listing  appears  in  Listing  10. 

Code:  Lines  6-12  Define  the  battlefield  and  dimension 

the  hill  arrays 

Lines  14-30  Read  the  hill  fitting  parameters  and 

convert  to  computing  parameters 

Lines  31-37  Reserve  and  read  in  the  ragged  array  LIST.H 

Lines  38-39  Determine  whether  there  are  any  forest  ellipses 

Lines  40-42  Dimension  the  forest  ellipse  arrays 

Lines  43-54  Read  the  ellipse  fitting  parameters  and 

convert  to  computing  parameters 
Lines  55-61  Reserve  and  read  in  the  ragged  array  LIST.C  . 
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VIII.  Some  Areas  for  Further  Development 

The  current  version  of  the  STAR  terrain  and  LOS  models  form  a working 
core  for  terrain  modelling  in  STAR. There  are  several  areas  in  which  work 
remains  to  be  done.  A preliminary  list  is  presented  here. 

1.  Incorporation  of  terrain  features  from  the  Arn\y  Mobility 
Model  patch  data  for  micro  terrain,  concealment,  vehicle 
limiting  speeds  (for  the  movement  model),  etc. 

2.  Development  of  a dynamic  smoke  model  in  conjunction  with 
detection  models  for  the  "dirty  battlefield". 

3.  Investigation  of  computer  assisted  terrain  fitting  to  speed 
the  hill  definition  task  and  possibly  result  in  better 
terrain  fitting. 

4.  Replacement  of  the  Iterative  Newton  algorithm  in  LOS  with 
a guaranteed  accuracy  approximation  using  splines  or  some 
other  appropriate  class  of  functions. 

Work  is  currently  underway  on  several  of  these  areas.  As  they  become  ready 
for  incorporation  into  the  main  STAR  model,  they  will  be  documented  separately. 
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